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We apply density functional theory to study the freezing of superfluid 4 He, charged bosons and 
charged fermions at zero temperature. We employ accurate Quantum Monte Carlo data for the 
linear response function in the uniform phase of these systems, a quantity that has different behav- 
ior for large values of the wavevector than previously assumed. We find that, as a result of this 
exact behavior, different approximations in the density functional theory of freezing that involve 
linear response, all fail to correctly describe the crystallization in three dimensions, while yielding 
satisfactory predictions in two dimensions. This demonstrates the shortcomings of the currently 
popular density functional approximate theories to describe 3d-freezing in the quantum regime. We 
also investigate the consequences of the exact asymptotic behavior of response functions on the form 
of effective interactions and polarization potentials in the electron gas, at small distances. 
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I. INTRODUCTION 

The modern density functional theory (DFT), which is 
employed in the theoretical investigations of freezing of 
both quantal and classical systems, is based on an exact 
correspondence between equilibrium one-particle densi- 
ties and external potentials.Ero In particular, if we denote 
by n(r) the one-particle density of the system (i.e. the 
statistical average of the one-particle density operator) 
the system can be characterized by an appropriate ther- 
modynamic potential which attains its minimum value 
for the correct (equilibrium) profile no(r). F° r the study 
of crystallization, the relevant thermodynamic potentials 
are the grand potential Q and the intrinsic Hclmholtz 
free energy F, the latter-, being a unique functional of 
the one-particle density.trH If (i is the chemical potential 
of the system at some temperature T and v ex t(v) is an 
arbitrary external potential, then the quantity: 



u] — F[n] — / rfrn(r)u(r) 



(1.1) 



where u(r) = fj, — v ext (r), is a minimum for given 
u(r) at the equilibrium density no(r). The quantity 
fl[u] = Q[no,u] is then the grand potential of the sys- 
tem. Clearly, the equilibrium condition reads as 



SF[r. 



5n(r) 



u(r). 



(1.2) 



«o(r) 



For vanishing external potential and fixed particle num- 
ber N, the intrinsic free energy F[n] is a minimum at 
the equilibrium density, with respect to variations of the 
shape of the density profile. It is customary to separate 
F[n] into a contribution from the noninteracting system 



under a suitable external potential that makes n(r) the 
equilibrium density, Fid[n], and an excess part F ex [n], i.e. 
F[n) = Fid[n] + F ex [n]. The determination of Fa is not 
complicated, for both classical and quantum systems: in 
the former case, Fid is known explicitly as a functional 
of the density!] In the latter case, the statistics appears 
explicitly in the construction of Fid, and for given ex- 
ternal potential one can construct both the equilibrium 
density and Fm in a straightforward manner. Therefore, 
the art of density functional theory amounts to the in- 
vention of approximate functionals for the excess part. 
In the classical regime, there has been extensive work 
in this direction during the last fifteen years .□ Relatively 
less has been done in the quantum regime, with which 
we are concerned in this work. 

The development of quantum DFT of freezing has fol- 
lowed two alternative routes: In one case,El suitable for 
finite temperatures, a mapping of the quantum particles 
into classical polymer rings is invoked; in the other, which 
is better suited for zem temperature, the Hohcnbcrg- 
Kohn-Sham formahsmoB is used, and the problem is, 
duced to a self-consistent band structure calculation! 
Here we follow the second approach, since we are inter- 
ested in T = freezing. In this case, F[n] is simply 
the intrinsic ground-state energy E[n\. The ideal part 
reduces to the kinetic energy of noninteracting 
particles To [n], and the remainder F ea; [n] is the excess 
energy E ex [n\. A brief summary of this formalism will 
be presented below. 

Within certain classes of approximate functionals, an 
essential ingredient for the practical implementation of 
this approach is the linear response function x( r ; n ) °f 
the fluid, or its Fourier transform x(q; n) where q is the 
wavevector magnitude and n is the average density. In 
particular, what is important is the 'quantum' direct cor- 
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relation function (dcf), i.e., the difference between the in- 
verse linear response functions of the interacting and non- 
interacting systems, ri) = x~ l {q\ n) — Xo n )- I n 
previous applicationaJTij it was assumed that this differ- 
ence is asymptotically vanishing (maybe in an oscillatory 
manner) for large YaJiWf s of the wavevector. However, re- 
cent exact results,00 and, associated Quantum Monte 
Carlo (QMC) calculations^^ Ej show that this is not the 
case: instead, the aforementioned difference approaches 
a negative constant as q — > oo. In this paper, we revisit 
the DFT of freezing, using the correct liquid-state input. 
We examine the perlarmance of the perturbative second- 
order theory (SOT)E£l and the nonperturbative modified 
weighted density approximation (MWDA) .113 For a vari- 
ety of systems, and irrespective of the range of the inter- 
action and the statistics (superfluid 4 He, charged bosons 
and fermions) we find that this exact large-g behavior 
has drastic consequences in three spatial dimensions: the 
crystal is predicted to be the stable phase for any den- 
sity. The SOT-functional is affected by this behavior 
most dramatically: it becomes unbounded from below 
as the density becomes more localized around the lat- 
tice sites and thus it has a minus- infinity minimum at 
a perfectly localized density. The MWDA, on the other 
hand, does not suffer from this extreme pathology: the 
MWDA-functional is bounded from below, but the (fi- 
nite) minimum of the energy always occurs for a mod- 
ulated (crystal) phase. In two dimensions, the effect is 
much less drastic, in the sense that for densities rele- 
vant to crystallization the SOT-functional continues to 
be bounded from below, yielding satisfactory predictions 
for the freezing of the electron gas. 

The rest of this paper is organized as follows: in Sec- 
tion II we present a summary of the DFT formalism; in 
Section III we survey the liquid-state input and discuss its 
implications on the behavior of the 'quantum' direct cor- 
relation functions, as well as on effective interactions — in 
the electrons gas; in Section IV we apply the SOT and in 
Section V the MWDA to the problem of freezing of differ- 
ent quantum liquids. Finally, in Section VI we summarize 
and conclude. 



II. QUANTUM DENSITY FUNCTIONAL 
THEORY OF FREEZING 

The quantum DFT formalism employed in this work 
has been presented in detail in Refs. [?||]. Here we give 
only an outline and refer the reader to the above papers 
for de tails . Wri ting E[n] = Xb[n] + E ex [n] and using 
Eqs. (1.1) and ( [1.2] ) we see that a necessary condition 



for equilibrium is 



5T [n] SE ex [n] 



8n(r) Sn(r) 



(Jt - v ext [r) 



(2.1) 



n (r) 



for the case of interacting particles. This is formally 
equivalent to the condition of equilibrium for noninter- 



acting particles (for which E[n] = Tq[ti\) under the influ- 
ence of an effective external potential 



Veff(r) 



5E R 



5n(r) 



Vext(r). 



(2.2) 



Note that the effective potential is itself a functional 
of the one-particle density, through the dependence of 
E ex [n] on n(r). Therefore, one is faced with a self- 
consistency calculation which in practice proceeds as 
follows: an initial guess is made for the density pro- 
file, which yields an initial form for the effective poten- 
tial. Then the one-particle Schrodinger equations (Kohn- 
Sham equations) 



- — V 2 + v eff (r) Vi(r) = e^(r) 
2m J 



(2.3) 



are solved, yielding the eigenfunctions ipi(r) and the as- 
sociated energy eigenvalues £j. From the former, a new 
one-particle density is constructed through 



n(r) =J2mMr)f 



(2.4) 



where n, are the occupation numbers suitable for the 
given statistics (Bose or Fermi). The new density serves 
for the construction of the new effective potential, and 
the cycle is continued until a self-consistent solution 
has been found. Once the self-consistent orbitals ipi(r) 
and the associated eigenvalues Si and density no(r) are 
known, the ideal kinetic energy To is given by 

T = ^n 4 |dr^(r)(-|lv 2 )^(r) 

i 

= y~] njEj - / drn (r)v eff (r). (2.5) 

i 

The formulation presented above is exact, provided v- 
representability holdsJIS Approximations enter through 
the excess energy functional S ra M which is not known 
exactly. In the following subsections we present two com- 
mon schemes which both rely on the knowledge of the 
second functional derivative of this functional with re- 
spect to the density at the uniform limit. This quantity 
is in turn directly related to the density-density linear 
response function. 

A. Second-order theory 

One usual approximation is the so-called second-order 
theory (SOT) or quadratic approximation. Here, one ex- 
pands functionally the unknown functional about a uni- 
form fluid of density n; , keeping terms up to second-order 
only. Explicitly, 



E ex [n] = E ex (m) + dr 



5E e . 



Sn(r) 



Sn(r) 



2 



1 

+ 2 



drdr' 



S 2 E ex [n] 



Sn(r)Sn(r') 



6n(r)Sn(r'), (2.6) 



with Sn(r) = n(r) — ni and E ex {ni) the excess intrinsic 
energy of the uniform liquid, a function of nj. Due to 
the translational and rotational invariance of the liqu id, 
the first functional derivative in the rhs of Eq. (2.6) is 
just a position-independent constant, equal to the excess 
chemical potential of the homogeneous liquid. The sec- 
ond functional derivative is a function of |r — r'| only; 
both depend on m, of course. We define, from now on, 



5 2 E ex [n] 



5n(r)Sn(r') 



-K(\r-r'\;rn). 



(2.7) 



The function K(r; n) is the excess part of the linear static 
inverse response function, of the homogeneous liquid, and 
can also be expressed aa3 



K(r;n)=x 1 (r; n) - x 1 (r; n), 



(2.8) 



where x _1 ( r ! n ) an d Xo ( r 'i n ) are the functional inverses 
of the density-density static linear-response functions of 
the interacting and noninteracting liquid, respectively. 

Such an approximation is not a priori guaranteed to 
have any validity, since there is no 'small parameter' guid- 
ing the expansion. Its widespread use is due on the one 
hand to practical limitations, as third- and higher-order 
functional derivatives of E ex [n] are poorly known even in 
the uniform phase, and on the other hand in the relative 
success that it has had, at least in the classical regime, 
in predicting the freezing parameters of simple liquids 
The function K(r; ni) is formally the quantum analog of 
the classical Ornstcin-Zcrnicke direct correlation function 

(dcf)B 

We set v ex t(r) = from now on. In the quadratic ap- 
proximation for the excess part of the energy functional, 
the effective potential which enters in the Kohn-Sham 
calculation is periodic with Fourier components 



^e//(Q) = Sn Q [ 



x 1 {Q;ni)] 



(2.9) 



where Q is a reciprocal lattice vector (RLV) of the given 
lattice and Shq is the Fourier component of the periodic 
function Sn(r) = n(r) — raj, and x _1 (q; ni) is the Fourier 
transform of the function x (i", ni). 

For systems of neutral particles, the choice of the den- 
sity ni of the reference liquid is arbitrary, although the 
usual choice is to consider a liquid at the same chemical 
potential as the solid. Moreover, for a Bose system at 
T = the kinetic energy of independent particles van- 
ishes in the uniform-limit. Thus, the difference between 
the grand potentiala of the solid and the liquid is:E£l 



AQ[n] = T [n] - - / / drdr' K{\r - r'\;ni)Sn{r)Sn{r') 

V ~ 



V 



\n Q \ 2 K(Q;ni). (2.10) 



Q#0 



In Eq. (2.1C), V is the volume of the system, n s is the 
average density of the solid, K(q; n) denotes the Fourier 
transform of K(r; n) at wavevector q, and uq is the 
Fourier component of the periodic density at RLV Q. 
In practice, one changes fi (or, equivalently, ni) and min- 
imizes Af2[n] with respect to n(r). Freezing occurs when 
min{Af2[n]} vanishes. For min{Af2[n]} > (< 0) the 
liquid (solid) is stable. 

For systems composed of particles carrying a charge e 
and interacting via the Coulomb potential v c (r) = e 2 /r, 
the presence of a uniform, rigid, neutralizing background 
of opposite charge guarantees the stability of the system. 
The presence of the background imposes the constraint 
that the freezing transition now takes place at constant 
density (isochoric freezing). The relevant thermodynamic 
potential is now the total energy E[n] \ the phase with the 
lowest E[n] is the thermodynamically stable one. It is 
customary for such systems to separate the excess energy 
into a Hartree contribution and an 'exchange-correlation' 
contribution, i.e. to write 



E ex [n] = — 



, , ,Sn(r)Sn(r') „ . , 
drdr' ^ ' , ' + E xc [n], (2.11) 



where Sn(r) — n(r) — n and n is the average density. If 
we now define 



S 2 E xc [n] 



Sn(r)5n(r') 



= -K xc (\r-r'\;m), 



then Eqs. (g/7j) and fl2.1l|) imply 
K{\r-r'\-m) = - Vc {\r-r'\)+K xc {\r- 



(2.12) 



;«,)• (2.13) 



In Fourier space, one writes the Fourier transform 
K xc (q; n) of K xc (r; n) as K xc (q; n) = v c {q)G(q; n), where 
V c (q) is the Fourier transform of the Coulomb poten- 
tial {v c (q) = 4ire 2 /q 2 in three dimensions and 2ire 2 /q 
in twOpdimcnsions) and G(q; n) is the so-called local field 
factorE3 Finally we have 



K(q;n)=v c (q)[l-G(q;n)}. 



(2.14) 



Due to the long-range nature of the Coulomb potential, 
the functional expansion of the energy of the inhomoge- 
neous phase can now be performed only about a liquid 
whose density n\ is equal to the average density n = n s 
of the solid. Using Eqs. (gj), ([D|), (p3Land ( gig ) we 
obtain the difference between the energyEj of the solid 
and the liquid phases as: 



AE[n] = T [n} - 



d + 2 



Ne F 



V 



2 ^ 

Q#0 



n Q \ 2 v c {Q)[l-G(Q;n s )}. (2.15) 
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Eq. ( [2.15| ) above is valid for fermions in d-dimensions 
with ep being the Fermi energy of noninteracting par- 
ticles in the liquid phase. For bosons, this equation re- 
mains valid with the omission of the second term in the 
rhs. 

As mentioned above, the lack of a small parameter in 
the functional expansion of the excess energy (at least as 
far as the freezing problem is concerned) has cast some 
doubt on the validity of the quadratic theory. This ob- 
servation has led to the development of a class of nonper- 
turbative approximations, which approximate the excess 
energy of the solid by that of a liquid. The density of 
the latter is a weighted average of the true density of the 
solid. Of particular interest, due to its computational 
simplicity and its success in describing bulk freezing for 
certain model systems in the classical regime, is the modi- 
fied weighted tLensity approximation (MWDA) of Denton 
and AshcroftjLZl presented in the following subsection. 



B. Modified weighted density approximation 

The MWDA amounts to the approximation of the ex- 
cess energy of the modulated system by that of a uniform 
system at a weighted densitytL3'BI with the latter being 
evaluated as a weighted average over the real density of 
the crystal in a self-consistent way. In other words, one 
writes 



E, 



Ef x WDA [n] 



Ne(fi), 



(2.16) 



where e(n) is the excess energy per particle of a uniform 
liquid of density n. The effective density n is evaluated 
as a weighted average over the spatially- varying density 
n(r) of the crystal and is defined by 



1 

N 



drdr' n(r)n(r')w(r — r'; h), 



(2.17) 



where the weight function w(r — r'; n), which depends on 
the weighted density itself is determined by requiring that 
the MWDA functional is exact to second-order in a func- 
tional expansion around a uniform liquid. The derivation 
of the expression for the .weight function has been pre- 
sented in detail elsewhereoEZI and so here we show only 
the final results which read as 



if n 

w{ " r]fl) = ~27(^) l^ r; ^ + v e " {fl - 



and 



n = n s + 



£ 

Q^O 



n Q \ 



-K(Q;n)] 
2e'(h) ' 



(2.18) 



(2.19) 



The effective potential for the MWDA is readily calcu- 
lated as! 



v eff(r) = e(n) + e'(n) 



Sn 
Sn(r) ' 



(2.20) 



and the corresponding expression in Fourier space, which 
is necessary for the solution of the MWDA-Kohn-Sham 
equations can be found in Ref. ^|. The MWDA excess en- 
ergy functional is exact to second order in a functional ex- 
pansion about a reference liquid, but also includes contri- 
butions from all higher orders. In this sense, the MWDA 
is a nonpcrturbative approximate scheme for the calcu- 
lation of the excess part of the energy. 



C. The Gaussian ansatz 

The self-consistent solution of the Kohn-Sham equa- 
tions is sometimes avoided by taking advantage of the 
fact that, in the solid phase, the particles are well local- 
ized around the lattice sites. This leads to the introduc- 
tion of the following Gaussian ansatz. One constructs 
normalized Bloch orbitals ipi,.(r) from a singleGaussian 
per site, 4>(r) = (2a/n) d / 4 e~ ar , according to:E3 



V>k(r) 



1 



s/NlW) 



R 



e ik-R e -a(r-R) 



(2.21) 



where {R} is the set of Bravais lattice vectors and 

P m (k) = J2R m e lWR - aR2/2 . (2.22) 



R 



After some algebra, we arrive at the following explicit 
expressions for the noninteracting kinetic energy To and 
the Fourier component of the density hq : 



and 



where 



Tain] 



'Q 



N— [da-a 2 n 2 ] 
2m 



M2 



-T 



P 2 (k) 



-Y 



P (k-Q/2) 
Po(k) ' 



(2.23) 



(2.24) 



(2.25) 



Eqs. fl2.23D-fl2.25D above, are valid for fermions; a de- 
notes the number of particles in each occupied orbital: 
(j = 1 for spin-polarized and a = 2 for unpolarized parti- 
cles. The k-sums extend over the occupied orbitals only. 
For bosons, we have to put all t he pa rticles in t he same 
orbital, k = 0. In this case Eqs. ( 2.23 ) and ( [2.24 ) remain 
valid with the identification: 



A*2 



P 2 (0) 

p (or 



Pq(Q/2) 
Po(0) 



(2.26) 



Substitutin g the appr opria te expression for To and tiq 
into Eqs. ( |2.10 ) or (2.15) above, one directly obtains 
the difference of the appropriate thermodynamic poten- 
tial between the solid and the liquid, within the SOT. In 
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t he M WDA, the additional self-consistent solution of Eg . 
( [2. 19 ) is required to get the excess energy of Eq. ( 2.16 ). 
In both cases one ends up with differences of thermody- 
namic potentials as a function of a, n s and n\. One then 
varies a (and n s for neutral particles) until a minimum 
is found. By repeating the procedure for different values 
of ni one can determine the phase diagram of the system 
at hand. 

We are going to present results obtained mainly 
through the use of the Gaussian ansatz, rather than the 
full self-consistent calculation. The reason is that, if a 
minimum exists when the Gaussian ansatz is employed, 
then the full calculation can only yield a lower minimum 
since the class of Gaussian densities is only a subclass 
of all the possible profiles. Since our calculations yield 
too low minima, the self-consistent calculation is for most 
purposes redundant. Moreover, the Gaussian ansatz al- 
lows for analytical estimates of th e m agnit udes of the 
ideal and excess terms in Eq. ( 2. lOf ) or ( 2.15 ). 

The excess liquid-state linear static inverse-response 
function K(r; n) plays, evidently, a central role in the 
implementation of the approximate schemes presented 
above. In the following Section we discuss the form and 
asymptotic behavior of this function. 



III. LIQUID-STATE INPUT, QUANTUM DIRECT 
CORRELATION FUNCTION, AND EFFECTIVE 
INTERACTIONS 

For classical liquids, the Fourier transform of the dcf 
is related to the experimentally measured structure fac- 
tor S(q) by a simple algebraic relationjij by virtue of the 
fluctuation-dissipation theorem. For quantum systems, 
on the other hand, the theorem relates dynamical quan- 
tities, and the relation between static quantities is not 
simple any more.0 As a result, various approximations 
for the static linear response function x(q) have been de- 
veloped. 

Superfluid 4 He is a test case. This is a fluid of neutral 
particles whose interactions can .be accurately described 
by the so-called Aziz potential.EJO In the absence of 
accurate data for x(q), one often resorts to the Feyn- 
man approximation to obtain a relation between x(q) 
and S(g)Ji3 which reads as 



XF{q;ni) = Xo{q;ni)S 2 (q), 



(3-1) 



where Xoil'^i) — —Amni/h 2 'q 2 is the static susceptibil- 
ity of the ideal boson gas. The ensuing approximate 
dcf Kp(q;ni) has been employed in density- functional 
theories-, of freezing of 4 He (Refs. fl0|]2"5| ) or Bose hard 
spheresu albeit with an appropriate 'rescaling' which was 
employed in an empirical way. This rescaling has been 
avoided in a recent density-functional study of quan- 
tum hard-sphere freezing.o However, accurate data for 
x(q) have now, been obtained from diffusion Monte Carlo 
calculations &3 
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FIG. 1. The function — niK(q;m) r(in mRy) of super- 
fluid 4 He as obtained from simulations ,lj for three different 
fluid densities. Solid line: rti = 0.02622 A -3 ; dashed line: 
= 0.02186 A" 3 : dash-dotted line: n, = 0.01964 A -3 . 



In Fig. ^ we show plots of this accurate direct corre- 
lation function for three different densities of the.liquid. 
A comparison with the Feynman approximation!!] shows 
immediately that whereas the latter has an oscillatory 
behavior about zero, the exact dcf is negative for almost 
all values of q > 2 A -1 . An additional important differ- 
ence concerns the large-q behavior of the dcf. Although 
the Monte Carlo data are limited to values q < 4 — 6 A -1 , 
exact theoretical calculations show that the q — > oo-limit 
of —K(q,ni) is a negative number, and not zero as the 
Feynman approximation implies.E3 In particulate, the re- 
sponse function x{Q'i n i) is given for large q by:EHI 



x(q;ni) = - 



Aran] 



H 2 q 2 



3H 2 q 2 



(KE)+0(q- i ) , (3.2) 



From Eqs. (2.8), (3.2) and using the result Xo{<i;ni) 



—Amni/h q z for the static susceptibility of the ideal bo- 
son gas, we obtain 



mK(oo;m) = — -(KE), 



(3.3) 



where (KE) is the expectation value of the kinetic energy 
in the liquid phase. These features of the exact dcf have 
important consequences on the performance of DFT's of 
freezing, as will be shown below. 

Charged fermions or bosons are another example of 
quantum liquids. The former is just the usual system 
of electrons in a uniform background (jellium) and the 
latter is a model system of spinless particles of electronic 
charge e and mass m in a background, but obeying Bose 
statistics. A natural length scale for these systems is the 
so-called Wigner-Seitz radius ro defined as the radius of 
a sphere which contains, on average, one particle, i.e. for 
a system of density n in d-dimensions we have 
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FIG. 2. The function ~niK(q;ni\ 
bosons as obtained from simulations: 
dashed line: r a = 50; dash-dotted line 
r s = 160. 



in mRy) of charged 
solid line: r s = 20; 
s = 100; dotted line: 



47rrj] 



(d = 3), 



2). 



(3.4) 



A convenient dimensionless measure of the density is 
r s = ro/ao, where ao is the Bohr radius. A widely used 
scheme to relate the local field factor G(q) with the struc- 
ture factor has been introduced by Singwi et a/E3 and is 
denoted by STLS. This has been employedJn DFT's of 
freezing of jellium in a number of cases.Ercl An impor- 
tant feature of the STLS scheme is that in the limit of 
large-g the local field factor G(q) approaches unity and 
this impli es tha t —K(q; m) approaches zero in that limit 
[see Eq. ( 2.1 4| ) ] . In this respect, the STLS scheme for 
systems of charged particles has the same features as the 
Feynman approximation. However, it has been shown 
axge-g limit, G(q) goes like q 2 in 



exactly that in thsplar 

three dimensions ;II2HI3 moreover it can .easily be shown 
that it scales like q in two dimensions. lj Inj-narticular, 
for charged bosons at d = 3 it is known that,E3 for large 

q, 



G(q; ni) 



2(KE)q 2 

iq{keY 



;(i-<?(o)) 



0(q- 



16((KE) 2 ) 
5h 2 u; 2 



(3.5) 



where lu p i — \J Annie 2 jm is the plasma frequency and 
g(0) is the value of the pair distribution funct ion o f the 
liquid g(r) at zero separation. From Eqs. (2.14) and 
(|3.5[) we find once more 



r H K(oc;n l ) = --(KE), 



(3.6) 



as in Eq. ( p. 3D above. 

For fermionSrin three dimensions, the large q local field 
factor reads astil 



G(q; ni) 



2{{KE) - (KE) )q 2 



+ 



16({(KEf) 



+ «(l-fl(0)) 



((KE) 2 ) ) 



5h 2 



16((KE) 2 - (KE)j 
9h 2 cu 2 



+ 0(q- 2 ), (3.7) 



where (• • -)o denotes a noninteracting average, and the 
coefficient of the q 2 term — the difference in the kinetic 
energy per particle between the interacting jaxid the non- 
interacting system — is a positive quantity.U Not e that 
the differences between Eq. (3.7) and Eq. ( |3.5[ ) arise 
from the different momentum distributions of the non- 
interacting Fermi and Bose systems. Using Eqs. ( p. 14 ), 
(|3.7|) we finally obtain, 



niK{oo; n{) 



\{{KE) - (KE) ) 



(3.8) 



In Fig. ^| we show the direct correlation function of 
charged bosons for a number of different densities as 
obtained from Quantum Monte Carlo simulations]i3 In 
Fig. H we show the same function for fully polarized 
charged fermions at r s = 100, which has also been ob- 
tained from QMCEj In both cases, it is clearly seen that 
at large values of q the function —K(q; n{) tends to a 
negative constant. In the system of point charged parti- 
cles, by virtue of the virial theorem, this constant may 
be expressed most simply as 



- niK(oo; n{) = - 



2 d{r s E) 



3 dr s 



(3.9) 



with E = e c (r s ) the correlation energy per particle, for 
fermions, and E = e(r s ) the energy per particle, for 
bosons. 
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FIG. 3. The function —niK(q;ni) (in 

mRy) of spiBppolarized charged fermions, as obtained from 
simulations,E3 at r s — 100. 
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FIG. 4. The function — mK(q;m) (in mRy) 
of spin-pola*ized charged fermions in 2d, as obtained from 
: 40. 



simulation E at r s 



In two dimensions, the situation is quite similar. For 
fermions, using the assapptotic behavior ot-the static lin- 
ear response function,El it has been showncS that the lo- 
cal field factor scales linearly with q, as q — > oo, namely 



Gfe „ = ffi£a^w, +1 _ 9(0)+ o ( o 



2ire 2 ni 



From Eqs. ( |2.14| ) and ( |3.10| ) we obtain 



n l K(oo;n l ) = -((KE)-(KE) ). 



(3.10) 



(3.11) 



In Fig. [| we show the direct correlation function 
of fully polarized electrons in 2 dimensions, near freez- 
ing, i.e., at r s — 4ft, as obtained from Quantum Monte 
Carlo simulations.E3 Again, the saturation of K(q, ni) to 
a constant which may be conveniently expressed as 

- mK(oo; m) = -((KE) - (KE) ) = (3.12) 

dr s 

with e c (r s ) the correlation energy per particle — is evi- 
dent. We note that the large q behavior of —niK(q,ni) 
for all the systems considered above, is given by 



niK(q,ni) = --((KE) - (KE) ) 
+0(q- d+1 )+0(q- 2 ), 



(3.13) 

^sterns 
X M that 



and evidently for the nonintcracting B 
(KE) q = 0. In fact one may easily she 
Eq. (3.13) above is valid for any quantum liquid inter- 
acting with pair potentials, both in in 3 and 2 dimensions, 
provided the second term on the rhs is only retained for 
Coulombic systems (1/r interaction) in 2 dimensions. 

The short-wavelength behavior of —K(q; nj) described 
above, implies that in real space the function —K(r; n{) 



has a delta-function contribution at the origin with neg- 
ative weight, as is clear from Eq. ( |3 .13 ). We shall there- 
fore define a regular def Kfi(q;ni), decaying to zero as 
q — > oo, by setting 

K(q;m) = K R (q;m) + -^((KE) - (KE) Q ). (3.14) 
nid 

This implies in real space 

K(r; m) = K R (r; m) + U^n^r)^ 1 

= K R (r;ni)+Ks(r;ni), (3.15) 
with the strength of the singular part Ks(r; ni) given by 



U (m) = -((KE) - (KE) Q ) > 0. 
d 



(3.16) 



Before investigating the consequences of this unexpected 
behavior of —K(r;ni) on the density functional theories 
of freezing, we shall pause here to briefly discuss its impli- 
cations on effective interparticle interactions in the liquid 
phase. As an example we shall consider spin unpolarized 
electrons in 3 dimensions. 

Within the dielectric formalism, the number and spin 
linear response functions of the normal electron fluid may 
be cast in a mean field, RPA-likp-form, by defining ap- 
propriate polarization potentials.^ In the static limit to 
which we shall restrict here, the number and spin re- 
sponse functions read respectively 



x(q) 



Xo(q) 



and 



Xs(q) 



i-v*( q )xo(q) 
xo(q) 



i-v a ( q )xo(q) 



(3.17) 



(3.18) 



with hb the Bohr magneton and V s (q) and V a (q) 
the symmetpc-and asymmetr ic p olariza tion potentials, 
respectively.E3l23 From Eqs. ( |2.S| ) and ( [2.14 ) it follows 
that 



V s (q) = 

with, G s (q,ri[) 
setEl 



K(q,m) =v c (q)[l-G s (q,m)], (3.19) 
= G(q,ni). In a similar fashion one can 



V a (q) = -v c (q)G a (q,n l ), 



(3.20) 



which defines the asymmetric local field factor G a (q, ni), 
whose behavior for large q is easily attained from the 
known asymptotic expansions of \s (<z)E3 as 

G a (q, m) = G s (q, m)-l + 2g(0) + 0(q- 2 ). (3.21) 

Interparticle polarization potentials for pairs of elec- 
trons with parallel or antiparallel spin projections are 
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readily obtained from their symmetric and asymmetric 
counterparts V s (q) and V a (q) viaE3 



V»%) = V s (q)±V a (q) 

= v c [l-G s (q,n l )TG a (q,n l )}, 



(3.22) 



where the upper sign corresponds to era ' = || an d the 
lower to aa' =tl- For large q, from Eqs. (3/7) and ( 3.21 ) 
one obtains 



and 



Vf°\q) = 2g(0)v c (q) + 0(q^). 



n 



(3.24) 



Eq. ( 3.23| ) implies the presence in (r) of a term 
-2U (ni)S(r)/ ni , with E7„(n,) > 0, given by Eq. pel) 
with <i = 3. On the other hand from Eq. ( 3.24[) one 
obtains that, for r -> 0, K p , ol (r) = 2.g(0)e 2 /r. This looks 



quite strange at first, as one would naively expect that at 
short distance effective interelectronic interactions should 
be essentially Coulombic. In fact, polarization potentials 
are not effective potentials, though at times this is not ap- 
preciated. We should also mention that in the approach 
of Refs. [32"|]35| the polarization potentials were assumed 
regular at the origin, V£°}(0) — e 2 q aa >, with the screen- 
ing wavevectors q aa i of the order of the Fermi wavevector 
q F . 

Effective electronic interactions have been defined for 
the elections gas by Kukkonen and Overhauser a long 
time ago,E-3 using the polarization potential method but 
keeping into account particle indistinguishability. Ac- 
cording to this study effective two-body electron-electron 
interactions may be written as 



V^{q) = v c (q) [1 + A acr ,(q)], 



(3.25) 



where 



A CTCT <(<z) 



v c (q){l-G s (q)] 2 Xo(q) 
l-V e {q)[l-G'(q)] X o(q) 
v c (q)[G a (q)] 2 Xo(q) 
l + v c (q)G"(q) X o(qy 



± 



(3.26) 



with the upper (lower) sign corresponding to parallel (an- 
tiparallcl) spins. The l arg e q beha vior of A CTCr / (q) is easily 
obtained from Eqs. r-(3.7), (3.21), and from the known 
iorEll of 



asymptotic behaviorLj of the Lindhard function 



Xo(q;ni)q^oo = - 



2 n2 



h q 



(3.27) 



One finds that A^(g) vanishes as q 2 for q — > oo, while 



27 



d(r s e c ) 



(3.28) 



with e c (r s ) the correlation energy per particle, in Ry- 
dbergs, of the electron gas. Thus V^(r) — e 2 /r for 
small r while the effective interaction between parallel 
spin is very slightly reduced with respect to the bare 
Coulomb repulsion, Vj-f(r) = r y(r s )e 2 /r, with j(r s ) — 
1 + Aff(oo) < 1. In particular, in the metallic regime 
one obtains_£rom the known equation of state of the elec- 
tron gasBEl j(r s ) = 0.99 and 0.98 for r s = 2 and 5. 
Thus, as we anticipated, effective interactions do remain 
essentially Coulombic at short distances. 



IV. SECOND-ORDER THEORY 

A. Three dimensions 

We begin with the application of the second-order the- 
ory (SOT) to the freezing of superfluid 4 He. Experimen- 
tal results on the system show that 4 He crystallizesEJ at 
a liquid density n; = 0.0260 A~ 3 into an hep-solid of den- 
sity n s = 0.0287 A~ 3 . We employ the Ga ussia n ansatz 
for a fee-crystal density and apply Eq. ( 2.1C ) for the 
evaluation of the grand-potential difference between the 
solid and the liquid. The value of K (g, n{) at q — which 
enters in this calculation is related to the energy per par- 
ticle e(ni) of the liquid via the 'compressibility sum rule', 
namely 



k{0-m) = 2e'{n l )+me"{n l ) 



(4.1) 



where the primes denote differentiation with respect to 
the argument. For the quantity e(ni) we use an analytic 
fit based on accurate diffusion Mont e Car lo data.Ea 

We try to minimize Af2[n] [Eq. (2.10)] with respect 
to a for a variety of different values of (n s , n{). As can 
be seen in Fig. || for the pair of values which are close to 
those for which freezing occurs in experiments, AO[n] has 
apparently no minimum] it keeps getting lower without 
bound as the localization increases. In the same figure 
it can be seen that for ni much lower than the freezing 
value, Af2 has a very negative local minimum at strong 
localization (large values of a), i.e. the solid is predicted 
to bee too stable. With reference to Fig. 0, note that at 
freezing one would expect aa 2 w 2, in order to obtain the 
correct 'quantum' Lindemann ratio 7 ~ 0.3. (7 is the ra- 
tio of the root mean square deviation about a site to the 
nearest neighbor distance.) On the contrary, the minima 
shown in the figure are at aa 2 ~ 10, implying a value 
of 7 oc l/\/a which is too small by about a factor 2, be- 
ing essentially classical. Unfortunately, the lack of Monte 
Carlo data for the def at large wavevectors does not al- 
low us to examine the limit of strong localizations, since 
as a grows we need more and more shells of RLVs into 
the sum of Eq. ( 2.10| ) in order to achieve convergence. 
Nevertheless, it is clear from the shape of —K(q; m) (see 
Fig. [I]) that the overestimation of the stability of the solid 
is brought about by the fact that —K(q, ni) is negative for 
all values of the RLVs; this way, the contribution from the 
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excess part of the energy, which is becoming lower with 
increasing localization, dominates over the contribution 
from the ideal energy, which grows with localization, to 
yield a total energy which is too low. We will make this 
statement more quantitative shortly. 



< 
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FIG. 5. Grand potential differenceEl [Eq. (2.10)] between 



a He fec-solid and a liquid at the same chemical potential, for 



different pairs (n s , ni 
n 3 = 0.0287 A" 3 , n, '■■ 
~ 3 0.0216 A" 



n, 



0.019 A -3 . Here 



in the second-order theory. Solid line: 
= 0.0262 A -3 ; dashed line: n s = 0.0287 
3 ; dash-dotted line: n s = 0.0275 A" 3 , 
, a = 2.556 A. 



Next, we look at the SOT-freezing of charged bosons, 
using the dcf-simulation results of Ref. 15. The system 
is known to undergo Wigner crystallization into a bec- 
solidt3 at r s — 160±10. Once more, we emplo y the Gaus- 
sian ansatz and try to minimize AE[n] [Eq. (2.15)] with 
respect to a at various different values of r s . Some of the 
results are shown in Fig. [|. The quantity AE[n] is clearly 
unbounded from below, i.e. the absolute minimum lies at 
infinite localization, where the value of A_E[n] is minus 
infinity. There is a local negative minimum at or 2 , « 3 
for r s = 50. This corresponds to the correct quantum 
Lindcmann ratio, and one might argue that the SOT of 
freezing may only make sense for modulations that are 
not too large (i.e., moderate values of a) and near the 
freezing density. Even so, the predicted freezing density 
would be overestimated by a factor of about 30. 

It becomes clear, therefore, that the SOT suffers from 
the pathology of producing an unbounded functional. 
Within the framework of the Gaussian approximation, 
this feature can be clearly understood as follows. Take a 
solid whose lattice constant is a and consider the strong- 
localization limit, i.e. a = aa 2 S> 1. In that case we 
have, with ex cellen t accuracy, \ii = and /io_= 1 for all 



Q's [see Eq. ( [2.26] )]. In e?-dimensions, Eq. (2.23) gives 
T (a 



(1+7?) 



h z d 

2ma 2 



a, 



(4.2) 



where n s — (1 + rf)ni (ry = for isochoric freezing) 



On 
the 



the other hand, the excess energy contribution 
sum o f the terms beyond Tq on the rhs of Eq. (2.10) or 
Eq. ( 2.15 )1 m &y be conveniently broken into two con- 
tributions originating respectively from Kn(q;ni) and 
Ks{q;ni). The first contribution is most easily treated 
in reciprocal space, where it takes the form 



Vm 



x ^ e -(Qa) 2 /4o n;/?i?(Q;n;); (4 3) 
Q#0 



and manifestly tends to a constant for large values of a. 
The second contribution is evaluated in real space, to 
leading order, as 



Vm 



5»1 



Uo(ni) 
2Vn 2 



dr{Sn(r)) 2 



2nia d ir d / 2 



(4.4) 



using the fact that, for a 3> 1, the density reduces to 
a superposition of nonoverlapping normalized Gaussians 
(2a/n) d / 2 exp{-2ar 2 }, one per site. The lack of a lower 
bound for d = 3 can be now easily understood: since the 
ideal energy scales like a and the excess like —c?l 2 for 
large a, it is then clear that their sum will be dominated 
by the — a 3 / 2 -term and will be unbounded from below. 
The analysis presented above is valid also for fermions. 
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FIG. 6. Ground-state energy differenceB [Eq. (2.15)] be- 
tween a charged boson bec-solid and the liquid, versus local- 
ization at three different average densities, as obtained from 
the second-order theory. Solid line: r s = 160; dashed line: 
r s — 100; dash-dotted line: r s — 50. 

If a minimization within the restricted space of Gaus- 
sian profiles fails to yield a finite minimum, then the 
global minimum of the unrestricted Kohn-Sham (KS) 







scheme, which cannot be higher, will also be minus infin- 
ity. This does not exclude, however, the possibility of ob- 
taining, by means of solving the KS-equations, some local 
minimum at a moderate value of the localization; in the 
KS-scheme there is no 'localization parameter' of course, 
but the Lindemann ratio, for example, can be used as 
a measure of the spatial extent of the one-particle den- 
sity around a lattice site. This possibility is particularly 
interesting because it could be argued that the perturba- 
tive character of the SOT immediately limits its validity 
to weakly-modulated density profiles. In order to pursue 
this line, "Wia-have also performed the full, self-consistent 
calculatiorJjj for both charged bosons and spin-polarized 
electrons, within the SOT, at the densities for which the 
dcf is available (see Figs. § and || above). However, no 
local minimum was found, in either case. The large-g 
behavior of the dcf and the number of space dimensions 
render the SOT pathological in d = 3. 



0.02 



and 




FIG. 7. Ground state energy difference between triangu- 
lar-solid and polarized liquid for electrons in 2d as function of 
the density. AE/N (in mRy) was calculated within the SOT. 
The solid squares are calculated points, with the line just a 
guide for the eye. 



AE ex (a;r s ) 



N 



1 2 d(r s e c ) _ /a y 

2 s dr. \ a J 



(4.6) 



where a is the lattice constant of the given crystal struc- 
ture and e c the correlation energy per particle in Ryd- 
bergs. Thus, at strong localizations 



AB(5;r s 



N 



1 2 d(r s e c )\ _ 
2 Ts dr. 



(4.7) 



For polarized fermions (a = 1) in 2d t he a vailable QMC 
datEO'cHl show that the coefficient in (4.7) is positive for 
r s < 59. Therefore, the SOT-functional remains bounded 
from below for values r s < 59. This is encouraging, given 
that the polarized electron gas in two dimensions crys- 
tallizes into a triangular lattice at a value of r s which is 
considerably smaller than this 'stability limit'. 



0.35 



0.3 




0.25 



FIG. 8. Lindemann ratio 7 around freezing in the trian- 
gular 2d Wigner crystal, as predicted by the SOT. The solid 
squares are calculated points, with the line just a guide for 
the eye. 



B. Two dimensions 



In two dimensions, Eqs. ( |4.2| ) an d (4.4) show that 
both the ideal and excess term scale as a at strong local- 
izations, the former with a positive and the latter with 
a negative coefficient. Therefore, the absolute values of 
the respective coefficients are crucial in determining the 
existence of a lower bound for the sum of the two terms. 
Exp res sing e nergi es in R ydbe rgs and making use of Eqs. 
(pd|, (O), Q3.16D , and fl3.12| ) we find for isochoric tran- 
sitions {rf = 0): 



T (a;r s 



N 



= 25 



(4.5) 



We have thus performed a full Kohn-Sham calcula- 
tion with the accurate liquid state input shown in Fig. 
^, using a-,Blane wave basis set as explained at length 
elsewhere. Qfl We have systematically checked conver- 
gence with respect to the plane wave cutoff and the num- 
ber of fc-points in the Brillouin zone. As it can be seen in 
Fig. § at r s =40, where we have the dcf from QMC, El the 
solid is still unstable, though its energy is only 6 micro- 
Rydbergs higher than that of the polarized liquid. On 
the ground that the explicit-dependence of G(q/qF',r s ) 
on r s should be very weak£j we have neglected it al- 
together to perform the calculations at the other val- 
ues of r s , using therefore the available local field factor 
G(q/ qF]T s = 40). This treatment predicts freezing from 
the polarized fluid at r s = 42 which agrees within error 
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bars with the QMC prediction of Tanatar and CeperleyE3 
r s = 37 ± 5 and is within two error bars from a more 
recent QMC predictioneS r s — 34 ± 4. We have also 
evaluated the Lindcmann ratio 7, which is shown in Fig. 
H as function of r s near freezing. We find 7 = 0.33 at 
r s = 42, to be compared with an accepted value for quan- 
tum freezing of about 0.3. We may thus conclude that 
the SOT is capable of predicting freezing in two dimen- 
sions with good accuracy. 

We have also investigated the effect of using a Gaus- 
sian ansatz for the Bloch orbitals (see, e.g., Sec. |II C| ), 
as opposed to the full Kohn-Sham calculation presented 
above. As expected, this yields solid energies that are 
slightly higher, predicting freezing at a lower density, i.e., 
at r s — 50 with 7 = 0.25. Though this is well out of the 
range predicted by QMC (see above) it is still a substait. 
tial improvement with respect to an earlier predictionE3 
of r s > 80, using the Gaussian ansatz and an approxi- 
mate local field factor. 

We should also mention results recently obtained using 
another DFT scheme,Ej in which a local density approxi- 
mation (LDA) to the total energy of the^modulated phase 
is augmented by gradient correctionsca (GCDFT). This 
approach yields crystallization at r s = 31, but with a 
density which appears to be very little modulated. In fact 
we have repeated such calculations, solving the equiva- 
lent one-orbital self-consistent Kohn-Sham equations in 
full. We reproduce the r s — 31 found in Ref. |d] and we 
find 7 = 0.369 which is very close to the uniform limit of 
7 = 0.373. To give a more direct idea of what this means 
in terms of localization, we may look at the minimum in 
the density profile at half distance between a site and one 
of its nearest neighbors, in units of the on-site density, 
5 = n(r nn /2)/n(Q). Here r nn is the nearest neighbor 
distance. We find that the GCDFT predicts S = 0.92 
at freezing, to compare with 6 = 0.30, which we have 
obtained within the SOT, and an exact value which is 
likelv,tp be even smaller. With respect to a conventional 
LDAEliJ in which the noninteracting kinetic energy is 
treated without approximation, the GCDFT is introduc- 
ing an overestimate of the kinetic cost of a modulation. 
For the small modulations predicted by the GCDFT, this 
may easily be checked by comparing, for instance at the 
first RLV of the triangular crystal, the exact 2d noninter- 
acting response functiorO'EEl with the one corresponding 
to the GCDFT, Xo°(q) = -(<7m/2nh 2 )/{l + q 2 /2q 2 F ). 
One might be tempted to argue that, for small modula- 



tions, the approach of Rot, 41 would be more consistent 



than a conventional LDAc.3 in that it treats all the com- 
ponents of the energy on the same footing. However, 
the quality of the resulting density profile, which is in- 
deed very poor, pointing to a weakly first-order if not a 
second-order transition, contradicts such a conclusion. 

Returning to the effects of the large-g behavior of the 
quantum dcf on freezing we may conclude that these are 
far less drastic in two dimensions than in three. Notice, 
however, that if one tried to apply the quadratic theory 
for systems with r s > 59, one would obtain, also in two 



dimensions, the erroneous answer that the stable phase 
is a crystal with infinite localization. This demonstrates 
that the quadratic theory gives a reasonable description 
of solids whose thermodynamic parameters are not far 
away from the freezing ones. Deeply inside the region of 
thermodynamic stability of the solid, the SOT loses its 
validity, even in those cases where it succeeds in predict- 
ing freezing. 

Having concluded the discussion of the quadratic the- 
ory in two and three dimensions, we now proceed with 
the nonperturbative approach, i.e. the MWDA. 



V. MODIFIED WEIGHTED DENSITY 
APPROXIMATION 

In this Section we will examine the behavior and per- 
formance of the MWDA for the case of systems of charged 
particles. The general analysis will show that regard- 
less of statistics, the MWDA yields a functional which 
is bounded from below, i.e. at the limit of large localiza- 
tions the energy difference between the solid and the liq- 
uid tends to +00 and not — 00 as in the case of the SOT. 
Then, we will present the application of the MWDA to 
the case of charged bosons, for which the availability 
of liquid-state input allows us to perform the MWDA- 
calculation. We still find the stability of the crystal to 
be overestimated, nevertheless. 

The presence of a singular term in the dcf K(r] n{) im- 
plies the existence of a similar t erm i n the weigh t function 
w(r;n), as is clear from Eqs. ( |2.1g| ) and ( 3.1 5| ) . With a 
straightforward analysis which closely parallels the one 
developed in the previous section for the excess energy 
one is led to the conclusion that for strong localizations 
(large values of the parameter a), the weighted density 
is given to leading order in a by 



U (n) 



-2he'(h)n d / 2 



J/2 



(5.1) 



suggesting that n grows with a as we shall demonstrate 
shortly, provided that —he'(n) > 0. We have verified that 
indeed —he'(h) ^^2Jdlf s e'(f s ) > for all the systems 
considered below.E3C3ll3 We shall examine the behavior 
of n in two and three space dimensions separately. 



A. Charged fermions in two dimensions 



On account of Eq. (5J) above, let us assume that n 
diverges with a and therefore that r s — the Wigner ra- 
dius in units of Bohr radii corresponding to the effective 
density h — goes to zero in the same limit. Using the def- 
initions of the previous section we may eliminate h in 
favor of f s to obtain 



1 



if 



:€c{fs)Y C . 

■a, 



f s e'(fs) r 2 



(5.2) 
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with the prime denoting differentiation with respect to f s 
and C a constant which depends on the structure chosen 
for the solid. Here e c and e are respectively the correla- 
tion and the excess energy (exchange plus correlation) of 
the 2 dimensional electron gas. ._. 
For small r s the excess energy is given byc3 



-Af\ 



B - Df, Inf. 



(5.3) 



with A, B, D positive constants which depend on the spin 
polarization. The dominant, — f^-term is the exchange 
energy and the remainder is the correlation energy e c (f s ). 
Using the above equation and keepi ng o nly the leading 
terms in the ratio appearing in Eq. (5.2) one obtains at 
once 



r s = {Ar 2 s /BCY' 6 a 



1/3^-1/3 



(5.4) 



Now, using Eq. ( 2.16 ) and si nce the excess energy of the 
liquid scales like — f" 1 , Eq. (5.4) yields 



E™ DA {r s ,&) 



N 



= -\E 2d (r s )\a 



~i/3 



(5.5) 



5»1 



Since the ideal term scales as a, it dominates over the 
excess one at large a and thus the total energy tends to 
plus infinity at the limit of strong localization. Notice 
that in the SOT the excess term scales as —a, whereas 
here only as —a 1 / 3 . The MWDA makes the dependence 
of the excess energy on the localization parameter a lot 
weaker than the SOT. We will see that this is also the 
case in three dimensions. 



B. Charged bosons and fermions in three dimensions 

In three dimensions, eliminating n in favor of f s in Eq. 
(P) yields 



j_ = ~{r s e c (f s ))' C -3/2 
f 3 f s e'(f s ) rf 



(5.6) 



where e c and e have the usual meaning for fermions, 
whereas for bosons e c coincides with e — the total energy 
per particle. Again, C is a constant depending on the 
structure assumed for the solid. As in two dimensions, 
we are led to assume that f s —> as a — > oo and there- 
fore we shall retain in this limit only leading terms in Eq. 

Charged bosons. As we have already mentioned 
e c (f s ) = e(f s ) in this case and as r, — * 0, we have to 



dominant order e(f s 
obtain 



-0.8031f: 



-3/4 



Ry 



Thus we 



(5.7) 



The energy e(f s ) now scales as — f s 3 ^ 4 , thus 



N 



= -\E b (r s )\a 



^3/8 



(5. 



The MWDA-excess energy per particle scales only as 
-d 3 / 8 as opposed to -o?l 2 in the SOT. Thus, the ideal 
energy which is linear in a dominates for strong localiza- 
tions, and the MWDA-functional is free of the pathology 
of the SOT, i.e. it does have a lower bound. 

The actual calculations that we carried out were for 
the case of bosons only; however, the same analysis can 
be carried out for fermions, and the results for this case 
are presented below. 

Charged fermions. The first few terms in the expansion 
of the excess grcoiJid-state energy of the electron fluid for 
small r s read asE3 



e(f s ) = -Bf' 1 +rinf s - A + 0(f s ), 



(5.9) 



where all constants are positive. Once more, the term 
proportional to — f^ 1 is the exchange energy and the re- 
mainder is the correlation energy. Using Eqs. (5.6) and 
([5.9|) we obtain as f s — > 



i ... M rci _ 3/2 

— =r s \\nr s \ — — c?' 2 . 
ri B r-i 



(5.10) 



To leading order as a — > oo we immediately obtain 



D{r s ) 



a- 3 / 8 

[lnd] 1 /^ 



(5.11) 



with D(r s ) = 



1/3TC} 1 / 4 . Finally, since the ex- 
the MWDA- 



cess energy per particle scales as —f s 
functional obeys the scaling 



E^^{r s ,a) 



N 



\Ef(r s 



dl 1 / 4 ^ 3 / 8 



(5.12) 



5»1 



Thus, the MWDA is free of the unboundedness problem 
also for fermions. It is interesting that the scaling of the 
exce ss ener gy is now dependent of the statistics [see Eqs. 
(|5.8|) and ( 5.12 ) above], though very weakly, due to the 
logarithmic dependence in a present for fermions. This 
is at variance with the prediction of the SOT, where the 
same scaling was found and it appears intriguing. In fact, 
naive considerations would suggest that the excess energy 
should scale in the same way for bosons and fermions 
at the strong localization limit since, in this case, each 
particle is confined to its own cell and statistics becomes 
unimportant. 

The existence of a lower bound for the MWDA- 
functional is an improvement over the behavior of the 
SOT. However, this property guarantees neither the exis- 
tence of a minimum at nonzero localization nor its correct 
location and behavior in terms of changes of the average 
density. If, for example, the total energy is monotonically 
increasing as a function of a, then the only minimum will 
occur for the uniform liquid. On the other hand, it is pos- 
sible that a minimum always exists, for any value of the 
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average density and is lower than the liquid one; in this 
second case, we are led to the erroneous prediction that 
the crystal is stable at all densities. In the following sub- 
section we show the results of the full MWDA-calculation 
for charged bosons and we find that, in fact, this second 
scenario materializes. 
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FIG. 9. Ground-state energy differencel-9 between the the 
charged boson bcc-solid and the liquid, versus localization at 
three different average densities, as obtained in the MWDA. 
(a) Solid line: r s = 20; dashed line: r s = 50. (b) r a = 100. 



C. MWDA-calculation for charged bosons 

We have im plem ented the MWDA-self-consistency 
condition [Eq. ( 2.1 9| ) ] for the case of charged bosons for 
which there exist sufficient simulation data for the local 
field factor for a range of densities varying from r s = 10 
to r s = 160 (Ref. |l^). We have used an analytic fit 
to the equation of state obtained from simulationEj We 
limit our study to the charged boson liquid because for 
polarized fermions the only available Monte Carlo data 
are for r s = 100 and the implementation of the MWDA 



requires the knowledge of the dcf of the liquid over a wide 
range of densities. 

We have carried out the MWDA-calculation with the 
Gaussian ansatz for three different values for the average 
density, namely r s — 20, 50 and 100. The results are 
shown in Fig. It can be seen immediately that, unlike 
the SOT, the MWDA gives minima of AE/N for finite 
values of the localization parameter a. Moreover, the 
trends of these minima are correct: they get deeper and 
also move towards stronger localization as the average 
density decreases. However, according to simulationsEj 
the liquid is stable for r s < 160± 10, whereas the MWDA 
gives a lower energy for the bcc-solid for values of r s as 
low as 20. Thus we can say that although the MWDA is 
already much better than the SOT, it still predicts a solid 
that is too stable. Even at high densities the MWDA- 
functional has a global minimum for a modulated phase. 



VI. CONCLUSIONS 

The implementation of the correct liquid-state input 
in a density-functional approach to the freezing of quan- 
tum liquids brings about a remarkable new result, namely 
that in three dimensions the standard SOT suffers from a 
lack of a lower bound and is always strictly minimized for 
infinitely localized solids. As we have already mentioned, 
one might argue that the perturbative character of the 
SOT limits its validity to weakly modulated density pro- 
files and thus local minima for finite localization should 
suffice. However, as we have demonstrated above for He 
and charged bosons, such minima — when they exist — still 
yield an incorrect description of freezing, predicting sta- 
bility of the solid well inside the region where the system 
in fact is liquid. Recourse to non-perturbative theories 
including certain classes of higher order terms, such as 
the MDWA, does not help much in practice. One gets 
rid of the extreme pathology of the SOT theory in that 
the resulting functional is bounded from below, which is 
certainly satisfactory. However even the MDWA predicts 
the crystal to be the stable phase deep into the region 
where the liquid should be stable. 

The asymptotic analysis in the localization parame- 
ter a that we have carried out above for the MDWA 
is easily generalized to the weighted density approxima- 
tion (WDA),e3Ij once it is realized that the weight func- 
tion w(q; n) has in this case the same large q limit as 
in MWDA, and therefore the same singular term in real 
space. One obtains the same scalings as discussed in Sec. 
[V] above, and therefore bounded functionals. Whether 
the predictions of the WDA for freezing will be any bet- 
ter than those of the MWDA remains to be investigated, 
though we 4ftubt it. We may mention that the conven- 
tional LDAtm also brings about bounded functionals, as 
one obtains almost at once. The scaling is the same as 
for the MWDA for electrons in 2 dimensions and charged 
bosons in 3 dimensions; it is different for electrons in 3 di- 
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a 1 / 2 . Again, the LDA is making the solid too stable 



mensions, for which the exchange-correlation energy goes 
like 

in three dimensions. u 

The situation in two dimensions appears specular to 
the one summarized above for three dimensions. In 
fact, earlier applications of the DFT theory of freezing 
with approximate liquid def gave good results, in three 
dimensions, Em while failing in two dimensions .E3 We have 
demonstrated above that the reverse is true, if more ac- 
curate liquid input is used, which obeys the exact large 
q behavior discussed in Sec. III. In particular, in two di- 



mensions the SOT provides bounded functionals, in the 
relevant region of density, and yields a good description 
on the freezing transition. 

The recently obtained accurate information on the 
liquid-state linear response functions gives rise, there- 
fore, to a new problem: our favorite approximate density- 
functional schemes which we have learned to trust from 
our experience on classical systems, seem to fail when 
applied to quantum systems in three dimensions. There 
is a need for reexamination of the current formalism of 
quantum density-functional theory of freezing and the 
development of approximate schemes which will be more 
appropriate to deal with the peculiarities of quantum 
systems. In this respect, the use of direct correlation 
functions possessing the correct asymptotic behavior is 
crucial, as is such a behavior that causes all the present 
troubles. At variance with the classical case we have seen 
that quantum functionals tend to predict excess energies 
(per particle) that negatively diverge at infinite localiza- 
tion. Though the MWDA produces a functional bounded 
from below, we speculate that the divergence of its excess 
part in this limit could still be incorrect, as the potential 
energy should remain finite, unless one can prove that it 
is the kinetic contribution to the excess energy which is 
bringing about this divergence. 
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